#include<vector>
#include<eigen3/Eigen/Dense>
#include<eigen3/Eigen/Sparse>
#include<cmath>
#include<algorithm>

class Function//函数类
{
public:
    virtual double operator()(double _x) = 0;//函数值
    virtual double diff(double _x){return 0;}//一阶导数
    virtual double diff_2(double _x){return 0;}//二阶导数
};

class Interpolation//插值类
{
public:
    virtual void solve() = 0;
    virtual double operator()(double _x){return 0;}//插值函数值
};

class Polynomial : public Function//多项式函数
{
public:
    Polynomial(std::vector<double> _coef, double _x0) : coef(_coef), x0(_x0){}
    double operator()(double _x)//f(x) = a0 + ... + an*(x - x0)^n
    {
        n = coef.size();
        double y = coef[0];
        for(int i = 1; i < n; i++)
        {
            double delta_y = coef[i];
            for(int j = 0; j < i; j++)
            {
                delta_y = delta_y*(_x - x0);
            }
            y = y + delta_y;
        }
        return y;
    }
private:
    int n;
    std::vector<double> coef;
    double x0;
};

class ppForm_interpolation : public Interpolation//分段多项式插值
{
public:
    ppForm_interpolation(Function &_f, std::vector<double> _Point, int _way, int _order) : f(_f), Point(_Point),
    way(_way), order(_order)//仅解决阶为1或3的情况
    {
        n = Point.size();
    }
    
    void solve()
    {   
        if(order == 1)
        {
            for(int i = 1; i < n; i++)
            {
                std::vector<double> ppf;
                double k;
                k = (f(Point[i]) - f(Point[i-1]))/(Point[i] - Point[i-1]);
                ppf.push_back(f(Point[i]) - k*Point[i]);
                ppf.push_back(k);
                C.push_back(ppf);
                Polynomial p(ppf, 0);
                P.push_back(p);//存储分段多项式
            }
        }
        if(order == 3)
        {
            Eigen::SparseMatrix<double> A(n,n);
            std::vector<Eigen::Triplet<double>> triplets;
            Eigen::MatrixXd b(n,1);
            b = Eigen::MatrixXd::Zero(n,1);
            for(int i = 1; i < n-1; i++)//存入对角线的大部分数据
            {
                double lamda = (Point[i+1]-Point[i])/(Point[i+1]-Point[i-1]);
                double miu = (Point[i]-Point[i-1])/(Point[i+1]-Point[i-1]);
                triplets.push_back(Eigen::Triplet<double> (i, i-1, lamda));
                triplets.push_back(Eigen::Triplet<double> (i, i, 2.0));
                triplets.push_back(Eigen::Triplet<double> (i, i+1, miu));
                b(i,0) = 3.0*miu*(f(Point[i+1])-f(Point[i]))/(Point[i+1]-Point[i]) 
                + 3.0*lamda*(f(Point[i])-f(Point[i-1]))/(Point[i]-Point[i-1]);
            }
            if(way == 1)//complete cubic spline
            {
                triplets.push_back(Eigen::Triplet<double> (0, 0, 1.0));
                triplets.push_back(Eigen::Triplet<double> (n-1, n-1, 1.0));
                b(0,0) = f.diff(Point[0]);
                b(n-1,0) = f.diff(Point[n-1]);
            }
            if(way == 2)//cubic spline with specified second derivatives at its end points
            {
                triplets.push_back(Eigen::Triplet<double> (0, 0, 4.0));
                triplets.push_back(Eigen::Triplet<double> (0, 1, 2.0));
                triplets.push_back(Eigen::Triplet<double> (n-1, n-2, 2.0));
                triplets.push_back(Eigen::Triplet<double> (n-1, n-1, 4.0));
                b(0,0) = 6.0*(f(Point[1])-f(Point[0]))/(Point[1]-Point[0])
                - f.diff_2(Point[0])*(Point[1]-Point[0]);
                b(n-1,0) = 6.0*(f(Point[n-1])-f(Point[n-2]))/(Point[n-1]-Point[n-2])
                + f.diff_2(Point[n-1])*(Point[n-1]-Point[n-2]);
            }
            if(way == 3)//natural cubic spline
            {
                triplets.push_back(Eigen::Triplet<double> (0, 0, 4.0));
                triplets.push_back(Eigen::Triplet<double> (0, 1, 2.0));
                triplets.push_back(Eigen::Triplet<double> (n-1, n-2, 2.0));
                triplets.push_back(Eigen::Triplet<double> (n-1, n-1, 4.0));
                b(0,0) = 6.0*(f(Point[1])-f(Point[0]))/(Point[1]-Point[0]);
                b(n-1,0) = 6.0*(f(Point[n-1])-f(Point[n-2]))/(Point[n-1]-Point[n-2]);
            }
            A.setFromTriplets(triplets.begin(), triplets.end());
            //A.makeCompressed();
            Eigen::SparseLU<Eigen::SparseMatrix<double> > lu;
            lu.compute(A);
            coef_x = lu.solve(b);
            for (int i = 1; i < n; i++)
		    {
                std::vector<double> ppf;
                double k;
                k = (f(Point[i])-f(Point[i-1]))/(Point[i]-Point[i-1]);
                ppf.push_back(f(Point[i-1]));
                ppf.push_back(coef_x(i-1));
                ppf.push_back((3.0*k-2.0*coef_x(i-1)-coef_x(i))/(Point[i]-Point[i-1]));
                ppf.push_back((coef_x(i-1)+coef_x(i)-2.0*k)/pow(Point[i]-Point[i-1],2));
                C.push_back(ppf);

                Polynomial p(ppf, Point[i-1]);
                P.push_back(p); 
		    }
        }
    }
    double operator()(double _x)//返回函数值
    {
        int k = 0;
        while(_x > Point[k])//找到x所在区间
        {
            k++;
        }
        if(_x == Point[k])
        {
            return f(Point[k]);
        }
        else
        {
            return P[k-1](_x);
        }
    }
private:
    std::vector<double> Point;
    Function &f;
    std::vector<Polynomial> P;
    std::vector<std::vector<double>> C;//分段多项式系数
    int n, way, order;
    Eigen::VectorXd coef_x;
};



class B_spline : public Function//B样条基函数
{
public:
    B_spline(int _order, std::vector<double> _t) : order(_order), t(_t){}
    double operator()(double _x)
    {
        std::vector<double> t1(t.begin(), t.end() - 1);//[t_i, t_i+1]
        std::vector<double> t2(t.begin() + 1, t.end()); //[t_i+1, t_i+2]
        B_spline B1(order - 1, t1);
        B_spline B2(order - 1, t2);
        if(order == 1)//1阶B样条
        {
            if(_x > t[0] && _x < t[1])
            {
                return (_x - t[0])/(t[1] - t[0]);
            }
            else if(_x >= t[1] && _x < t[2])
            {
                return (t[2] - _x)/(t[2] - t[1]);
            }
            else
            {
                return 0;
            }            
        }
        else//根据递归求阶不为0时的B样条
        {
            return(_x - t[0])*B1(_x)/(t[order] - t[0]) 
            + (t[order+1] - _x)*B2(_x)/(t[order+1] - t[1]);
        }
    }
private:
    std::vector<double> t;
    int order;//B样条阶数
};

class B_spline_interpolation : public Interpolation//仅支持线性样条和三次样条
{
public:
    B_spline_interpolation(Function &_f, std::vector<double> _Point, int _way, int _order) : f(_f), 
Point(_Point), way(_way), order(_order)//f为待插值函数，Point为插值点(不妨取等距的), way为插值方式, order为样条阶数
    {
        n = Point.size();
        double d = Point[1] - Point[0];//插值点距离
        if(order == 1)//插入额外点
        {
            double x0 = Point[0];
            double xn_1 = Point[n-1];
            Point.insert(Point.begin(), x0-d);
            Point.push_back(xn_1+d);
        }
        if(order == 2)
        {
            double x0 = Point[0];
            double xn_1 = Point[n-1];
            Point.insert(Point.begin(),x0-d);
            Point.insert(Point.begin(),x0-2*d);
            Point.push_back(xn_1+d);
            Point.push_back(xn_1+2*d);
        }
        if(order == 3)
        {
            double x0 = Point[0];
            double xn_1 = Point[n-1];
            Point.insert(Point.begin(),x0-d);
            Point.insert(Point.begin(),x0-2*d);
            Point.insert(Point.begin(),x0-3*d);
            Point.push_back(xn_1+d);
            Point.push_back(xn_1+2*d);
            Point.push_back(xn_1+3*d);
            // point.insert(point.begin(),x0-d);
            // point.insert(point.begin(),x0-2*d);
            // point.push_back(xn_1+d);
            // point.push_back(xn_1+2*d);
            // for(int j = 0; j < n+1; j++)
            // {
            //     std::vector<double> i(point.begin()+j, point.begin()+j+4);
            //     B_spline b(2, i);
            //     B2.push_back(b);//B-2样条
            // }
        }
        for(int j = 2-order; j < n+1; j++)
        {
            std::vector<double> i(Point.begin()+j-2+order, Point.begin()+j+2*order);
            B_spline b(order, i);
            B.push_back(b);//存储所需的所有B样条
        }   
    }

    void solve()
    {
        if(order == 1)
        {
            Eigen::SparseMatrix<double> A(n,n);
            std::vector<Eigen::Triplet<double> > triplets;
            Eigen::MatrixXd b(n,1);
            b = Eigen::MatrixXd::Zero(n,1);
            for(int i = 0; i < n; i++)
            {
                triplets.push_back(Eigen::Triplet<double> (i, i, 1.0));//令A的对角线为1
                b(i,0) = f(Point[i+1]);
            }
            A.setFromTriplets(triplets.begin(), triplets.end());//初始化系数矩阵
            A.makeCompressed();
            Eigen::SparseLU<Eigen::SparseMatrix<double> > lu;//LU分解
            lu.compute(A);
            coef_x = lu.solve(b);//计算Ax=b
        }
        if(order == 2)
        {
            double left = f(-4.5);
            double right = f(4.5);
            Eigen::SparseMatrix<double> A(n,n);
            std::vector<Eigen::Triplet<double>> triplets;
            Eigen::MatrixXd b(n,1);
            b = Eigen::MatrixXd::Zero(n,1);
            for(int i = 1; i < n-1; i++)
            {
                triplets.push_back(Eigen::Triplet<double> (i , i-1, 1.0));
                triplets.push_back(Eigen::Triplet<double> (i, i, 6.0));
                triplets.push_back(Eigen::Triplet<double> (i, i+1, 1.0));
                b(i,0) = 8.0*f(Point[i+2]);
            }
            triplets.push_back(Eigen::Triplet<double> (0 , 0, 5.0));
            triplets.push_back(Eigen::Triplet<double> (0 , 1, 1.0));
            triplets.push_back(Eigen::Triplet<double> (n-1 , n-2, 1.0));
            triplets.push_back(Eigen::Triplet<double> (n-1 , n-1, 5.0));
            b(0,0) = 8.0*f(Point[2])-2.0*left;
            b(n-1,0)=8.0*f(Point[n+2])-2.0*right;
            A.setFromTriplets(triplets.begin(), triplets.end());//初始化系数矩阵
            A.makeCompressed();
            Eigen::SparseLU<Eigen::SparseMatrix<double> > lu;//LU分解
            lu.compute(A);
            coef_x.resize(n+2);
            coef_2 = lu.solve(b);//计算Ax=b
            coef_x(0) = 2.0*left-coef_2(0);
            coef_x(n+1) = 2.0*right-coef_2(n-1);
            coef_x.segment(1,n) = coef_2;
        }
        if(order == 3)
        {
            Eigen::SparseMatrix<double> A(n+2,n+2);
            std::vector<Eigen::Triplet<double>> triplets;
            Eigen::MatrixXd b(n+2,1);
            b = Eigen::MatrixXd::Zero(n+2,1);
            for(int i = 1; i < n+1; i++)
            {
                triplets.push_back(Eigen::Triplet<double> (i, i-1, B[i-1](Point[i+2])));//A的三对角线
                triplets.push_back(Eigen::Triplet<double> (i, i, B[i](Point[i+2])));
                triplets.push_back(Eigen::Triplet<double> (i, i+1, B[i+1](Point[i+2])));
                b(i,0) = f(Point[i+2]);
            }
            if(way == 1)//complete cubic spline
            {
                std::vector<double> v0(Point.begin()+1, Point.begin()+5);
                std::vector<double> v1(Point.begin()+2, Point.begin()+6);
                std::vector<double> vn_1(Point.end()-6, Point.end()-2);
                std::vector<double> vn(Point.end()-5, Point.end()-1);
                B_spline b0(2, v0), b1(2, v1), bn_1(2, vn_1), bn(2, vn);
                triplets.push_back(Eigen::Triplet<double>(0, 0, (-3.0*b0(Point[3]))/(Point[4]-Point[1])));
                triplets.push_back(Eigen::Triplet<double>(0, 1, (3.0*b0(Point[3]))/(Point[4]-Point[1]) - (3.0*b1(Point[3]))/(Point[5]-Point[2])));
                triplets.push_back(Eigen::Triplet<double>(0, 2, (3.0*b1(Point[3]))/(Point[5]-Point[2])));
                triplets.push_back(Eigen::Triplet<double>(n+1, n-1, (-3.0*bn_1(Point[n+2]))/(Point[n+3]-Point[n])));
                triplets.push_back(Eigen::Triplet<double>(n+1, n, (3.0*bn_1(Point[n+2]))/(Point[n+3]-Point[n]) - (3.0*bn(Point[n+2]))/(Point[n+4]-Point[n+1])));
                triplets.push_back(Eigen::Triplet<double>(n+1, n+1, (3.0*bn(Point[n+2]))/(Point[n+4]-Point[n+1])));

                // triplets.push_back(Eigen::Triplet<double> (0, 0, (-3.0*B2[0](Point[3]))/(Point[4]-Point[1])));
                // triplets.push_back(Eigen::Triplet<double> (0, 1, (3.0*B2[0](Point[3]))/(Point[4]-Point[1])-
                // (3.0*B2[1](Point[3]))/(Point[5]-Point[2])));
                // triplets.push_back(Eigen::Triplet<double> (0, 2, (3.0*B2[1](Point[3]))/(Point[5]-Point[2])));
                // triplets.push_back(Eigen::Triplet<double> (n+1, n-1, (-3.0*B2[n-1](Point[n+2]))/(Point[n+3]-Point[n])));
                // triplets.push_back(Eigen::Triplet<double> (n+1, n, (3.0*B2[n-1](Point[n+2]))/(Point[n+3]-Point[n])-
                // (3.0*B2[n](Point[n+2]))/(Point[n+4]-Point[n+1])));
                // triplets.push_back(Eigen::Triplet<double> (n+1, n+1, (3.0*B2[n](Point[n+2]))/(Point[n+4]-Point[n+1])));
                b(0,0) = f.diff(Point[3]);
                b(n+1,0) = f.diff(Point[n+2]);//边界条件
            }
            if(way == 2)//cubic spline with specified second derivatives at its end points
            {
                triplets.push_back(Eigen::Triplet<double> (0, 0, 6.0/((Point[4]-Point[1])*(Point[4]-Point[2]))));
                triplets.push_back(Eigen::Triplet<double> (0, 1, -6.0/((Point[4]-Point[1])*(Point[4]-Point[2]))-
                6.0/((Point[5]-Point[2])*(Point[4]-Point[2]))));
                triplets.push_back(Eigen::Triplet<double> (0, 2, 6.0/((Point[5]-Point[2])*(Point[4]-Point[2]))));
                triplets.push_back(Eigen::Triplet<double> (n+1, n-1, 6.0/((Point[n+3]-Point[n])*(Point[n+3]-Point[n+1]))));
                triplets.push_back(Eigen::Triplet<double> (n+1, n, -6.0/((Point[n+3]-Point[n])*(Point[n+3]-Point[n+1]))-
                6.0/((Point[n+4]-Point[n+1])*(Point[n+3]-Point[n+1]))));
                triplets.push_back(Eigen::Triplet<double> (n+1, n+1, 6.0/((Point[n+4]-Point[n+1])*(Point[n+3]-Point[n+1]))));
                b(0,0) = f.diff_2(Point[3]);
                b(n+1,0) = f.diff_2(Point[n+2]);
            }
            if(way == 3)//natural cubic spline
            {
                triplets.push_back(Eigen::Triplet<double> (0, 0, 6.0/((Point[4]-Point[1])*(Point[4]-Point[2]))));
                triplets.push_back(Eigen::Triplet<double> (0, 1, -6.0/((Point[4]-Point[1])*(Point[4]-Point[2]))-
                6.0/((Point[5]-Point[2])*(Point[4]-Point[2]))));
                triplets.push_back(Eigen::Triplet<double> (0, 2, 6.0/((Point[5]-Point[2])*(Point[4]-Point[2]))));
                triplets.push_back(Eigen::Triplet<double> (n+1, n-1, 6.0/((Point[n+3]-Point[n])*(Point[n+3]-Point[n+1]))));
                triplets.push_back(Eigen::Triplet<double> (n+1, n, -6.0/((Point[n+3]-Point[n])*(Point[n+3]-Point[n+1]))-
                6.0/((Point[n+4]-Point[n+1])*(Point[n+3]-Point[n+1]))));
                triplets.push_back(Eigen::Triplet<double> (n+1, n+1, 6.0/((Point[n+4]-Point[n+1])*(Point[n+3]-Point[n+1]))));
                b(0,0) = 0.0;
                b(n+1,0) = 0.0;
            }
            A.setFromTriplets(triplets.begin(), triplets.end());//初始化系数矩阵
            A.makeCompressed();
            Eigen::SparseLU<Eigen::SparseMatrix<double> > lu;//LU分解
            lu.compute(A);
            coef_x = lu.solve(b);//计算Ax=b
        }
    }

    double operator()(double _x)
    {
        if(order == 1)
        {
            if(_x == Point[1])
            {
                return f(Point[1]);
            }
            int k = 1;
            while(_x > Point[k+1])
            {
                k++;
            }
            return coef_x(k-1)*B[k-1](_x) + coef_x(k)*B[k](_x);
        }
        if(order == 2)
        {
            if(_x == Point[2])
            {
                return f(Point[2]);
            }
            int k = 1;
            while(_x > Point[k+2])
            {
                k++;
            }
            // double B2_1 = 0;
            // double B2_2 = 0;
            // double B2_3 = 0;
            // B2_1 = (Point[k+2]-0.5-_x)*(Point[k+2]-0.5-_x)/(((Point[k+2]-0.5)-(Point[k]-0.5))*((Point[k+2]-0.5)-(Point[k+1]-0.5)));
            // B2_2 = (_x-(Point[k]-0.5))*((Point[k+2]-0.5)-_x)/(((Point[k+2]-0.5)-(Point[k]-0.5))*((Point[k+2]-0.5)-(Point[k+1]-0.5)))
            // +(((Point[k+2]-0.5)-_x)*(_x-(Point[k+1]-0.5)))/(((Point[k+3]-0.5)-(Point[k+1]-0.5))*((Point[k+2]-0.5)-(Point[k+1]-0.5)));
            // B2_3 = (_x-(Point[k+1]-0.5))*(_x-(Point[k+1]-0.5))/(((Point[k+3]-0.5)-(Point[k+1]-0.5))*((Point[k+2]-0.5)-(Point[k+1]-0.5)));
            return coef_x(k-1)*B[k-1](_x) + coef_x(k)*B[k](_x) + coef_x(k+1)*B[k+1](_x);
        } 
        if(order == 3)
        {
            if(_x == Point[3])
            {
                return f(Point[3]);
            }
            int k = 1;
            while(_x > Point[k+3])//找到_x所在区间
            {
                k++;
            }
            return coef_x(k-1)*B[k-1](_x) + coef_x(k)*B[k](_x) + coef_x(k+1)*B[k+1](_x) + coef_x(k+2)*B[k+2](_x);
        }
        return 0;
    }
private:
    std::vector<double> Point;
    std::vector<double> point;
    std::vector<B_spline> B2;//way=1时，矩阵需要用
    std::vector<B_spline> B;
    Function &f;
    int n, way, order;
    Eigen::VectorXd coef_x;//方程系数
    Eigen::VectorXd coef_2;
};



